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ABSTRACT 

We present new methods for lensing reconstruction from CMB temperature fluctu- 
ations which have smaller mean-field and reconstruction noise bias corrections than 
current lensing estimators, with minimal loss of signal-to-noise. These biases are usu- 
ally corrected using Monte Carlo simulations, and to the extent that these simulations 
do not perfectly mimic the underlying sky there are uncertainties in the bias correc- 
tions. The bias-hardened estimators which we present can have reduced sensitivity to 
such uncertainties, and provide a desirable cross-check on standard results. To test 
our approach, we also show the results of lensing reconstruction from simulated tem- 
perature maps given on 10 x 10 deg 2 , and confirm that our approach works well to 
reduce biases for a typical masked map in which 70 square masks each having 10' on 
a side exist, covering 2% of the simulated map, which is similar to the masks used in 
the current SPT lensing analysis. 

Key words: gravitational lensing - cosmic microwave background - cosmology: 
observations. 



1 INTRODUCTION 

Ongoing, upcoming and next-generation CMB experi- 
ments are able to measure arcminute-scale temperature 
anisotropics, which are perturbed significantly by gravita- 
tional lensing. Recently, several studies have reported the de- 
tection of lensing signals by reconstructing lensing fields in- 
volved in the CMB anisotropics, using th e cross-correlations 
between CMB and large-scale structure (ISrnith et al]|2007l ; 
iHirata et alj|2008l : iBleem et alJliail; ISherwin et al.ll2012t) 
or C MB maps alone ( Das et al.l l201ll ; Ivan Engelen et al.l 
l2012h . The measurement of the lensing effects through 
upcoming and next-generation CMB experiments will be 
a powerful probe of the properties of dark energy and 
massive neutrinos (e g.. iHul |2002l; iLesgourgues fc Pastorl 



direction n, axe given by 

9(n) = 8(n + d(h)) = Q(n) + d(n) • V9(n) + C(|rf| 2 ) . 

(1) 

The vector, d(n), is the deflection angle, and, in terms 
of parity, we can decompose it into two terms, these are 
the gradient (even pa rity) and curl (odd parity) modes 
|Namikawa et al.ll2012l ): 



d(n) = V<t>{h) + (*V)ro-(n) . 



(2) 



120061 : Ide Putter et al.|[2009l;lNamikawa et al.ll20ld). primor 



dial n on-Gaussianity (e.g., iJeong et al.||2009l; iTakeuchi et al.l 
20121). and cosmic strings (e.g.lYamauchi et al.l l201ll ; 



Namikawa et al.ll2012l ; lYamauchi et alJl2012T ) 



The distortion effect of lensing on the primary temper- 
ature anisotropics is expressed by a remapping. Denoting 
the primary temperature anisotropics at position n on the 
last scattering surface, O(n), the lensed anisotropies in a 
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where the symbol, *, denotes an operation which rotates 
the angle of a two-dimensional vector counterclockwise by 
90 degrees. Hereafter, we refer to <f> and zn as the scalar and 
pseudo-scalar lensing potential, respectively. 

Estimators for the lensing defle ction field have been 
derived by several autho r s (e.g . . IZaldarriaga fc Seliakl 



19991: ISefiak fc Zaldarriagal 1 19991; Irki 



Qkamotd |200~ 

Okamoto fc Hull2003l: IHirata fc Seliakl l2003bl : iBucher et al.l 
Namikawa et al. I I2OI2I ). These estimators all utilize 



2012; 



the fact that a fixed lensing potential introduces statistical 
anisotropy into the observed CMB, in the form of a correla- 
tion between the CMB temperature and its gradient. With a 
large number of observed CMB modes, this correlation may 
be used to form estimates of the lensing potential. The power 
spectrum of the lensing potential may in turn be studied by 
taking the power spectrum of these lensing estimates. 

When performing lens reconstruction on a realistic 
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dataset, there are two sources of bias which must be cor- 
rected for: 

(i) Statistical anisotropy from non-lensing sources such 
as a sky mask, inhomogeneous map noise, and asymmetry 
of the instrumental beam can be misinterpreted as lensing, 
generating a spurious lensing "mean-field". 

(ii) In the power spectrum of the lensing estimates, there 
can be a large contribution from the reconstruction noise 
which must be subtracted. 

With perfect statistical understanding of the CMB, fore- 
grounds, and instrumental response to the sky these sources 
of bias are always computable (although potentially only 
through the use of Monte Carlo simulations), and may be 
subtracted from the lensing potential and power spectrum 
estimates. The bias terms are often quite large in com- 
parison to the lensing signals of interest, however, the re- 
sult of which is that uncertainties in the power spectrum 
of the CMB fluctuations, the shape of the instrumental 
beam and its transfer function, the contribution from un- 
resolved (and therefore unmasked) point sources, and the 
instrumental noise level can all lead to problematic uncer- 
tainties for these bias terms. Multiple approaches have been 
proposed in the literature to mitigate these problems. For 
mean-field biases, most studies have focused on the issue 
of mask ing. One approach is to explicitly avoid m asked 
regions (|Hirata et al.ll2008l : ICarvalho fc Te'renoll201ll ). An- 
other is inpainting, in which masked regions are filled with 
simulated signal, th erefore reducing the large gradients at 
the m ask boundary (|Perotto et al.ll2010l ; iPlaszczvnski et all 
120121 ). Published analyses from the ACT and SPT experi - 
ments have used either source sub traction (|Das et al.l 201 lh 
or inpainting via Wiener filtering (jvan Engelen et aI.ll2012T ) 
to deal with resolved point sources, and apodization to re- 
duce spurious gradients at the survey boundary. The reduc- 
tion of reconstruction noise bias has also been discussed in 
the literature. For full-sky coverage, with homogeneous map 
noise and a symmetric beam the reconstruction noise bias 
may be estima t ed directly from the po wer spectrum of the 
map (|Hull200ll ; Ibvorkin fc Smithll2009l ). This approach has 
the added benefit of suppressing t erms in the covarian ce of 
the reconstructed power spectrum (|Hanson et al.ll2~01ll ) . For 
several specific choices of apodization/inpainting, it has been 
found that the fu ll-sky equations can be quite accurate eve n 
for cut-sky data (|Das et alj|201ll ; IPlaszczvnski et alj|2012t ). 
In principle, the reconstruction noise bias may be avoided 
entirely by taking the cross-spe ctrum of two estimators with 
independent noise realizations (|Hull200ll ). This is a more dif- 
ficult proposition than for power spectrum, where the only 
source of noise is instrumental and independent surveys of 
the same region of sky may usually be obtained. In the case 
of lensing, a large fraction of the reconstruction noise comes 
from the CMB fluctuations themselves, and so the construc- 
tion of estimators with independent noise realizations re- 
quires slicing the Fourier pla ne into d isjoint regions, such 
as the odd/even pari ty split (Hu 2001) or the in/out split 
(|Sherwin fc DastolOT ). There is usually a substantial loss of 
signal-to-noise associated with such splits. Additionally, for 
a realistic observation, some mode mixing is induced by the 
sky mask and instrumental response, and it is necessary to 
introduce buffer regions between the disjoint pieces of the 



Fourier plane from which the lensing estimates are formed, 
leading to further degradation of signal-to-noise. 

In this paper, we discuss new methods for constructing 
lensing estimators which have significantly reduced mean- 
field and reconstruction noise biases, with minimal loss of 
signal-to-noise. These "bias-hardened" estimators can be 
constructed in conjunction with any of the data processing 
methods (such as inpainting, apodization, inverse-variance 
filtering) described above. They are useful to deal not only 
with the complications induced by masking, but also with 
other effects such as noise inhomogeneity, beam asymmetry, 
and uncertainty in the primary CMB power spectrum. 

This paper is organized as follows. In Sec. [2] we briefly 
summarize quadratic estimators for lensing potentials from 
the CMB temperature, and the biases which they must be 
corrected for. In Sec. [3] we propose methods for avoiding 
these biases. In Sec. [4] we perform numerical simulations to 
test how well these methods work for several choices of map 
filtering and characteristic sky cuts. Section [5] is devoted to 
summary. 



2 LENSING RECONSTRUCTION 

In this section, we briefly review the general formalism for 
quadratic temperature lensing estimators on the flat-sky 
jHu fc Okamotdl2002l ; ICoorav et ai1l2005l ; iNamikawa et all 
l2012t ). 



2.1 Lensing potential estimators 

The statistical anisotropy introduced by lensing of the CMB 
is given by 

<e,_ i e i ) C MB = s t cf e + fir.**, ( 3 ) 

Here, to distinguish from the ensemble average over realiza- 
tions of the observed map including noise (■•■), we denote 
(• • -)cmb as the ensemble average with a fixed scalar and 
pseudo-scalar mode of lensing potential. The weight func- 
tion is given by 

flz. = Cf e £ Q X L + Cf t %t & x (£—£), (4) 

where x = <j> or zu, and the operator, Q x , is defined, for 
arbitrary two vectors, a and b, as 

aQ<j,b = a - b, a ro b = (*a) • b . (5) 

We note that the weight given in Eq. @ uses the lensed 
power spectrum instead of unlensed power spectrum, fol- 
lowing |LewisiFal] (|201ll ). 

This motivates the following generic form for a 
quadratic lensing estimator: 

Xi = xf - (xf) , (6) 

where the first term is the "standard" quadratic estimator 



(-) 



and a correction for the mean-field bias is given by the sec- 
ond term. Here Af is a normalization^, and Ql are the 
Fourier modes of a filtered sky map. 

1 In principle the normalization here should be a matrix, however 
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The filtered sky map O t, may be obtained in a variety of 
different ways. Inverse- variance (or "C~ ") filtering can be 
shown to minimize the reconstruction noise llHirata fc Seliakl 
l2003al ; ISmith et all 120071 ; lHanson fc Lewis! I2009I ). For full- 
sky coverage and homogeneous map noise, with a delta func- 
tion instrumental beam so that the data model is 

9(n) = e(n)+n(n), (8) 

where O is a (lensed or unlensed) CMB realization and n(n) 
is a noise realization, the corresponding diagonal filter is 
given by 

®L = -^e L (diagonal) . (9) 

The quantity Cf e , is the theoretical ensemble average angu- 
lar power spectrum of the observed sky, including the contri- 
bution from instrumental noise. This diagonal filter is simple 
to implement, and can be used even on the masked sky, with 
the penalty of large spurious gradients at the mask boundary 
|Hirata et al.ll200ct ). Intermediate between full C _1 filtering 
and the diagonal approximation is the approach of using an 
apodized sky mask to reduce the creation of spurious gradi- 
ents at the mask boundary. In this paper we will study all 
three approaches, in conjunction with bias-hardened estima- 
tors. In all of our discussion, we will assume that whichever 
filter is chosen, in regions far from any mask boundary it 
asymptotes to the form given in Eq. (J9]). 

On the full sky, with diagonal filtering, the estimator 
normalization may be determined analytically and is given 
by 

At -yw^j ■ (10) 

In more general situations this expression does not necessar- 
ily hold and the normalization must always be determined 
using Monte Carlo. It is still incredibly useful, however, if the 
effective normalization is close to that given by Eq. (|10[) . as 
this equation may be used to propagate uncertainties on the 
weight function of Eq. Q (for example, due to uncertainties 
in the instrumental beam transfer function) to uncertainties 
on the normalization without the need for expensive simu- 
lations. 

A numerical approach to fast computation of the generic 
quadratic estimator is to rewrite Eq. Q as a convolution of 
two maps, ©£, and on = iiCf e Q t (|Hu fc Okamotoll2002h : 

xf =A$f ■j0j s e L \i® x at-i.]. (11) 

With the convolution theorem, in real space this can be writ- 
ten as 

xl=A? J d 2 he- z£ - A e(h)[£Q x a(h)}, (12) 

where the quantities, Q(n) and a(n), are the real-space 
counterpart of 6v and ot£, respectively. Eq. (|12p means that 
the estimator can be computed by Fast Fourier Transform. 



for most realistic situations this is impractical and instead an 
effective normalization like the one above is used. 



2.2 Power spectrum estimators 

The power spectrum of the lensing potential may be stud- 
ied through the power spectra of the quadratic estimators 
above. This quantity probes the 4-point function of the 
lensed CMB, and can be usefully broken into disconnected 
and connected parts as 

<|i,T> = <|i?! 2 > D + <|i, c ! 2 >c. (is) 

The disconnected part, (■•■)d, contains the contributions 
which would be expected if Ql were a Gaussian random 
variable, while the connected part, (• • -)c, contains the non- 
Gaussian contributions which are a distinctive signature of 
lensing. The disconnected part represents the reconstruction 
noise bias discussed earlier, which must be accurately sub- 
tracted from Eq. (|13[) to obtain a clean measurement of the 
lensing signal. 

The disconnected bias is given by 

(M )n = -(A e ) J— J — 

x fe,i,fe,L'C'L,e-L,'C'£-L,L'> (14) 

where Cl,l' — {QlOl') is the covariance matrix of the fil- 
tered map. Given a model for this covariance matrix and 
a method to simulate Gaussian realizations of it, this dis- 
connected bias may be evaluated by Monte Carlo. If the 
covariance matrix is diagonal then it is tractable to eval- 
uate the disconnected bias analytically. For the filtering of 
Eq. ((9}, the disconnected bias is equal to the normalization, 
with A\ = {\x^\ 2 ) D = N*- {0 \ where A? is given by Eq. (fTO)). 

The power spectrum of the connected part of the 
quad ratic estimator is given on the full sky by l|Kesden et al.l 
|2003( ) 

<|x?| a >o = Cr + JVr (1) - (15) 

Here C% x is the potential power spectrum which it was our 
intention to reconstruct, while N^'^ is a nuisance term 
coming from the "sec ondary" lensing contractions of the 
trispectrum (|Hull200l| ). 



3 BIAS-HARDENED ESTIMATORS 

In this section, we propose methods to mitigate the mean- 
field and reconstruction noise biases discussed in the previ- 
ous section. 

3.1 Bias- reduced lensing estimator 

We begin by discussing a method to reduce the lensing 
mean-field. Our approach is straightforward; we construct 
a new estimator which is optimized to detect the source of 
the mean-field bias, and use this estimate to correct the lens- 
ing estimator accordingly. This allows the construction of a 
new hybrid lensing estimator which has intrinsically much 
smaller mean-fields than the standard one. 

To illustrate this approach, it is useful to consider the 
mean-field bias introduced by a specific form of statistical 
anisotropy. Here we will focus on masking, although the bi- 
ases from b eam asymmetry and inho mogeneous map noise 
are similar (|Hanson et al. I l2009l . l2010h . To start, consider a 
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modulation of the observed temperature fluctuations by a 
function W(n): 



9 mod (n) = W(h)e>(n). 



(16) 



From Eq. (|16[) . the multipole coefficients of observed temper- 
ature anisotropies are related to the underlying fluctuations 

a.s 

er d = J d 2 ne u - A e mod (n) = J (17) 

where we use the Fourier counterpart of the window func- 
tion: 



We = / d 2 he itfl W{h). 



(18) 



Introducing a mask function Me = 5e — We, we can rewrite 
this as 



er" = e. - / JU-m,_,.q,, 



(19) 



(2tt) 2 

The covariance matrix of the masked sky is then given by 

/AmodAmod\ c A©@ , \ 1 rx _ 

(OL Oj-i)CMB,n = OlC L + 2_, U,L x t 



+ f&Me + 0{M 2 ) , 



(20) 



where we define f t L — — (7® 
(• • -)cMB,n means the ensemble average over the realizations 
of CMB and noise with a fixed realizations of lensing poten- 
tials. At first order in M, with diagonal filtering, Eq. (|20[) 
leads to the following bias for the scalar lensing potential 
estimator: 



WcMB,n = 4>l + ■ 



d 2 L ft L f, 



LJt,L 



(2tt) 2 2Cf @ Cf, 



Me . (21) 



We see that masking introduces a mean-field which directly 
traces the mask Mi. Unlike the scalar lensing potential, the 
estimator for pseudo-scalar lensing potential is unmodified 
at first order in M. Masking does not introduce a large 
mean-field into the pseudo-scalar lensing potential. 

An estimator for the mask field, Me, can be constructed 
analogous to that for lensing: 



Mf = Uf 



(22) 



where A\ l is the same as Eq. (|10[) but using the mask weight 
function, f™ L . 

Now we consider the joint estimation of both the mask 
and lensing fields simultaneously. The standard quadratic 
estimator is biased by masking as Eq. (|21[) . Correspondingly, 
the naive mask estimator of Eq. 1|22|) is biased by lensing as 



<M|)cMB,n = Me + 



A 



d 2 L 



rM r4> 
Ji,L.Tt 



l.L 



(2tt) 2 2C? e Ce e , 



fa. (23) 



In matrix form, we can write (temporarily ignoring the 
mean-field corrections for both estimators) 



(WCMB.11 

(Mf) 

CMB,n 



l 



tie 



1 



Mi 



(24) 



where the ensemble average is taken over CMB and noise 



realizations and we define the response function R^' , for 

a,b = <j>,M, 

d 2 L fe,L,fe.L 



R a >" = Ai 



O) 2 2cf e qf e i| 



(25) 



By inverting Eq. p4[) , we obtain 



91 
Me 



tin ti,, 



1 



-ti p 



£f )cMB,n 



(Mf) 



CMB,n 

(26) 

A bias-free estimator for the scalar lensing potential may 
therefore be formed as 

(27) 



RY Me 



1 - tt e K e 

This lensing estimator has no mean-field contribution from 
the mask (up to the M 2 term in Eq. I20[) . Of course, in 
a practical situation the estimator of Eq. ([270 will not be 
completely free of mask bias. The M 2 term may produce 
a mean-field contribution, and also if non-diagonal filtering 
of the map is utilized then the response terms of Eq. ()25|) 
are only approximate. Even in such situations, however, the 
estimator of Eq. (|27[) should have a smaller mask mean-field 
than the standard estimator. We therefore refer to this as 
a "bias-reduced (lensing) estimator". In Sec. U we will test 
the behavior of this estimator for several choices of filtering 
and mask. The procedure above can be easily generalized to 
mitigate multiple sources of mean-field simultaneously. 

In the case of masking (as well as beam asymmetries 
and noise inhomogeneity) , the shape of the mean-field Mi 
is already known perfectly and we would also subtract the 
mean-field bias using the naive subtraction of Eq.©. How- 
ever, the bias-reduced estimator still has importance since 
it is the optimal unbiased estimator in the presence of mask 
field, i.e., the above estimator is uniquely determined by 
imposing the unbiased and optimal conditions on both the 
lensing and mask estimators, in the presence of both lensing 
and masking effects 0. 

Our intention for the construction of bias-reduced es- 
timators is to mitigate the effect of possible errors in our 
understanding of the mean-field. Even if we know the mean- 
field Me perfectly, there are other sources of possible error 
for the mean-field subtraction of Eq.((S}. Suppose, for exam- 
ple, that our analysis is performed using a slightly incorrect 
estimate Cf e (denoted by a cursive C) for the ensemble av- 
erage map power spectrum Cf e , with 

Sz, . (28) 



abb _ ago 

L'L — 



This propagates directly to an error ff L — — E_l — £|£-z, 
in the weight function of Eq. (|20j> . If the mean-field bias is 
removed by averaging over masked CMB realizations with 
power spectrum given by Eq. (|28|) . there will be an uncor- 
rected mean-field contribution given by 



TV'^Mi , (29) 
where !Zg b is the same as Eq. (|25[) but using the incorrect 



2 This is the same analogy of iNamikawa et al. but now 

we consider the lensing and mask fields, which do not separately 
estimate each other. Thus, we have to consider the response in- 
duced by the other, as shown in Eq . d251 l . 
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power spectrum. For the bias-reduced estimator, however, 
the uncorrected contribution is given by 



-M t . 



(30) 



It is in principle possible that the residual mean-field in this 
case is worse than for the standard approach, for example 
if is zero but TZ M ^ is not. If a source of non-lensing 

statistical anisotropy has sufficient resemblance to lensing 
to produce a significant mean-field, such a situation seems 
somewhat pathological however and we do not believe it to 
be common. A more likely situation is that of a calibration 
error, for example ff L = bfi' L , for some small coefficient 
b. In this case the residual mean-field will be completely 
avoided by the bias-reduced estimator. If bounds on E can 
be obtained for a specific experiment, the usefulness of the 
bias-reduced estimators can be explored using Eq. (|29|l and 
Eq. ()30p . In any case, agreement between standard and bias- 
reduced estimators provides a useful consistency test. 

3.2 Noise bias estimator 

We turn now to the issue of reconstruction noise bias, given 
by Eq. (|14|l . which is dependent on the covariance matrix 
of the filtered CMB modes Cl i,i. Similar to the case of 
the bias-reduced lensing estimator above, we suppose that 
we are in possession of an imperfect model C^i/ for the 
ensemble-average covariance matrix of the filtered CMB 
modes 

Cl,l* = Cl,l' ~ ^l,l> , (31) 

where T. L L i is an error matrix. An estimate of the recon- 
struction noise made by substituting C for C in Eq. (|14l) will 
have 0(E) contributions from the error matrix. We would 
therefore like to construct an estimator to determine the 
reconstruction noise bias more directly from the data. 

For full-sky coverage and diagonal filtering, where the 
covariance matrix is given by — 8 £ _ L iCl, this can 

be done simply, by replacing the ensemble average Cl in 
Eq. (|14[) with the ( realization-dependent) power spe ctrum 
of the filtered map l|Hull200ll : [Dvorkin fc Smithll2009h . This 
method of correcting the disconnected bias has the added 
advantage that it removes the largest off-diagonal contri- 
butions to the covariance ma trix of the power spectrum 
estimates (|Hanson et all [20 111 ). In more realistic situations 
where the covariance matrix has off-diagonal elements this 
procedure is not guaranteed to work, although for some spe- 
cific forms of filtering it has been found adequate ( Das et ah! 
l201ll : Ivan Engelen et alj|2012l ; IPlaszczvnski et alj|2012i ). 

Here we motivate a new approach which utilizes both 
data and the imperfect covariance. It is more robust than 
relying entirely on Cj, ; £/ , and also does not depend on the ac- 
curacy of full-sky equations which neglect any off-diagonal 
correlations due to masking, inhomogeneity of the instru- 
mental noise, etc. The method is again straightforward, we 
simply estimate the reconstruction noise bias using Eq. (|14p . 
substituting the imperfect model for one of the covariance 
matrices and the data itself for the other: 



<|i,|% = (A?) 2 ± 



d 2 L 



d 2 L' 



pX pX 

(2tt) 2 J (2n) iJt ' Llt ' V 



(32) 



This approach to removal of the disconnected bias emerges 
naturally when deriving optimal trispectrum estimators 
fro m an Edgeworth ex pansion of the CMB likelihood (see 
e.g. iRegan et all |2010l ) for a thorough derivation) . The cal- 
culation of the bias in this manner is only sensitive to uncer- 
tainties in the CMB covariance at 0(E 2 ), an improvement 
over an entirely model-based determination of the recon- 
struction noise. It also maintains the property of suppressing 
off-diagonal contributions to the covariance matrix of the re- 
constructed power spectrum. We refer to the approach above 
as the "noise bias estimator". 



4 NUMERICAL TESTS 

In this section, we test the usefulness of the bias-reduced 
estimator for several different choices of filtering and char- 
acteristic masking. As discussed in Sec. 13.11 we expect sig- 
nificant reduction of the mask mean-field bias using these 
estimators, which reduces the ability for some uncertain- 
ties in the primary CMB and instrument properties to leave 
residual biases in the estimated lensing potential. We do not 
test the noise bias estimator, as the only potential limitation 
of its usefulness is in the size of the E 2 matrix in Eq. (|31[) . 
and this must be evaluated on an experiment-specific basis. 
We will, however, compare the noise bias and normaliza- 
tion for the bias-reduced estimators to the approximation of 
Eq. pop as a cross-check on how closely they agree with the 
full-sky expectation for various choices of filtering, with or 
without the use of bias-reduced reduced estimators. 

For lensing reconstruction, we use 20 realizations of 
Gaussian unlensed and lensed temperature fluctuations. 
These fluctuations are simulated on a 10 x 10 deg 2 map, 
and the lensed maps are generated using the ray-tracing 
simulations described in Appendix [X] We generate simu- 
lated source masks by cutting iV m randomly located square 
regions with angular size r m on each side. In our analysis, 
we choose N m — 70 and r m = 10' or 20'. Note that the case 
with Nm — 70 and r m = 10' roughly corresponds t o the case 
of SPT lensing analysis (jvan Engelen et all 120121) . To con- 
sider experiments with high-angular resolution such as Po- 
larBear, ACTPol and SPTPol, and also to avoid contamina- 
tion by SZ and unresolved point sources, we assume a delta 
function instrumental beam, but truncate the temperature 
multipoles at £ max = 3000. We assume homogeneous map 
noise, with a level of 1 fj,K arcmin. The angular power spec- 
trum of the lensed/unlensed temperatu re and scalar lens - 
ing potential is computed using CAMB (|Lewis et al.ll2000l ). 
Throughout this section, our fiducial values of cosmologi- 
cal parameters are c onsistent with WMAP 7-year results 
(|Komatsu et all 1201 ll ); fl h = 0.046, fl m = 0.27, Ov = 0.73, 
h = 0.70, a 8 = 0.81, n s = 0.97. 



4.1 Filtering 

We use three approaches to filtering our simulated skymaps: 
the straightforward diagonal filtering of Eq. <(9j , on a masked 
map, diagonal filtering on the map with an apodized mask, 
and C _1 filtering. The apodization and C" 1 procedures are 
described in more detail below. 
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Figure 1. Real part of the Fourier counterpart of the apodizing 
function given in Eq. I|34|l with so = 0.0 (red solid line) compared 
with that of the top-hat function (so = 1.0; green dashed line). 
We choose a = 10 deg. 



4-1.1 Apodization 

One approach to reduce mode coupling from sky cuts, ofte n 
used in power spectrum estimation (e.g. jDas fc Bodell200ot ). 
is apodization; to smooth the mask somewhat so that its 
Fourier counterpart more closely resembles a delta function. 
To apodize the survey boundary, for example, we can use a 
window function given by (as in Eq. I16|) 



W s (x, y; s ) = w a (x; s )w s (y; s ) . 



(33) 



apodizing window function to the observed map, given as 

iv m 

W(h; s , to) 



Wi(s ,*o) 



W s (n;s )l[(l-W^(n-t j), 

i=l 

(36) 

The function, W s (n; so), is used to apodize the edges of the 
survey region, while the functions 1 — Wr^\(n;to), apodize 
the point sources. The factor Wi is given by 

r iv m 
W 1 (s ,t ) = / d 2 he lM W\h; so) ]](1-W§(n;t )), 

(37) 

and the functions, (h; to) [i = 1,2, ... , N m ), are defined 



W^(x,y,to) = w^(x - Xi;to)w^(y - yv,to) 



(38) 



with 



7T fr(l+to)- |fj 

2 bt Q 



\t\<b 

b<\t\< 6(1 + to) 
6(1 +t ) < \t\ 



(39) 

The positions (xi,yt) denote the position of i-th source 
mask, and 6 = r m /2. The parameter, to, indicates the 
size of the apodization region for each source mask, and, 
similar to the case of W B (n;so), the Fourier counterpart, 
W$(£;to) = / d 2 ne ufl W(^(n;to), has a sharp peak at 
£ ~ for large values of to- If both functions, W s (n; so) and 
Wf^\(n;to), are sharply peaked at £ ~ in Fourier space, 
the Fourier transform of W^n; so, to) can be approximated 
as a delta function. 



We will use a sine apodization function given by 



w B (s; so) — — x < 

Wx 



7T 1 

sm,-- 




■S() 



|s| < as 

aso < \s\ < a . (34) 
a < \s\ 



The parameter, so, indicates the width of the region where 
the apodization is applied, and the prefactor, Wi = 2a[so + 
2(1 — so)/7r], is used so that ds w s (s, so) = 1. In Fig. [1] 
we show the Fourier counterpart of the above function: 



w B (£;s ) 



j its s / \ 

as e w (s; so) 



(35) 



The Fourier counterpart with so = 0.0 has a high-contrast 
peak at £ = relative to that of the top-hat function 
(so = 1.0). This implies that the function W s (n;so) given 
by Eq. (|34[) with so = 0.0 would be a better choice to reduce 
mode coupling, compared to that with so = 1.0. From here 
forward, the parameter, so in all of our survey boundary 
apodization is set to so = 0.0. 

We now construct an apodization function for both the 
survey boundaries and detected point sources. Let us con- 
sider an observed map given on an S — [— a : a] x [— a : a] 
plane with Nm detected point sources, each of which we 
would like to mask. For simplicity, we will use a square mask 
function, with a length of r m on each size. We apply an 



4-1.2 C" 1 filtering 

The minimum-variance filtering which emerges from 
likelihood-based derivations of lensing estimators is known 
as C~ l filtering. For the data model of Eq. ((Sj the inverse- 
variance filtered multipoles, , are obtained by solving 



[l + C^iV^C 172 ] (C 1/2 0) = C 1/2 N- 1 Q , 



(40) 



where © is a vector whose components are O^, C is the 
covariance of the lensed or unlensed CMB anisotropics with 

^ee 



{C} e 



= S, 



(41) 



and N = (n^n) is the covariance matrix for the instrumen- 
tal noise. The noise covariance matrix in Fourier space is 
obtained from that in real space as 



AT 



(42) 



(43) 



where the pointing matrix, Y, is defined by 

{Y} fli .i ] = exp(m< • £j) . 

The mask is incorporated by setting the noise level of 
masked pixels to infinity, and therefore the inverse of the 
noise covariance in real space N to zero for masked pixels. 
The inversion of the matrix on the left-hand side of Eq. ((40} 
can be numerically costly, but may be evaluate d using con- 
jugat e descent with careful preconditioning (|Smith et alj 
l2007t) . 
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Figure 2. The power spectrum of the mask mean-field, using 20 Gaussian unlensed map realizations with the standard quadratic 
estimator (left) and the bias-reduced estimator (right). In each panel, we show the case with diagonal filtering (S-diag, BR-diag), 
apodization (S-apo, BR-apo), and C-inverse filtered map (S-inv, BR-inv). The black line shows the Monte-Carlo noise, 

N p(°)/ 20 . The 

multipoles are used up to £ m ax = 3000. The number and size of masks, N m and r m , are fixed with 70 and 10', respectively; the total 
fraction of masked area is ~ 2%. 
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Figure 3. Same as Fig[2] but for r m = 20'; the total fraction of masked area is ~ 8%. 



For our C _1 results, after applying the C" 1 filter we 
additionally apply an apodizing function to account for the 
survey boundary, given by Eq. (p4")l . 



4.2 Mean-field power spectrum 

We now proceed to our numerical results. We start by look- 
ing at the power spectrum of the mean-field for both the 
standard and bias-reduced estimators. 

In Fig. [2l we plot the power spectrum 



W 2 J 2tt 



1 

20 



(44) 



where xf'' is reconstructed from i-th realization of an un- 
lensed Gaussian map without mean-field subtraction, and 
Wi is derived (in analogy to Appendix [B| as 



W 2 



d hW (n; so, to) ■ 



(45) 



We construct source masks with N m = 70 and r m = 10'. To 
show the usefulness of bias-reduced estimator, we compute 

• the standard quadratic estimator (Eq. [6)), or 

• the bias-reduced estimator (Eq. 1270 . 

where the filtering is taken to be 

• the diagonal filter with no apodization of source holes 
(to = 0, denoted as S-diag) or 

• the diagonal filter with apodized source holes (non-zero 
to, denoted as S-apo), or 

• the inverse- variance filtered map (denoted as S-inv). 

As noted above, in all cases we use so = 0.0. The results for 
the bias-reduced estimator are prefixed by "BR" . 

It is clear that the mean field bias from the standard 
quadratic estimator is large particularly on large scales, 
500. For the standard quadratic estimator, the mean-field 
on large scales still has a large amplitude even with source 
apodization or C" 1 filtering. When we use the inverse- 
variance filtering, the mean-field contribution from source 
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Figure 4. The angular power spectrum of lensing estimator with Gaussian simulation, using the standard quadratic estimator with the 
subtraction of the mean- field bias (Eq. |7J] upper panel) and the bias-reduced estimator fEo. 1271 lower panel). The mean and error bars 
of angular power spectrum arc computed by 20 realizations of the simulation. 
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Figure 5. The angular power spectrum of the bias-reduced 
estimator for pseudo-scalar lensing potential, in the case with 
N m = 70 and r m = 10'. 



holes is suppressed significantly. This is because the at term 
in the estimator (c.f. Eq. Ill)) in the case of CT 1 filtering 
corresponds to Wiener filtering, which is able to reconstruct 
the component of the masked signal which is due to modes 
larger than the holes themselves. Most of the power in the 
CMB gradient is due to scales greater than r m = 10', and so 
this inpainting aspect of the C _1 filter significantly reduces 
the generation of spurious gradients near the source bound- 
ary. Even for C _1 filtering, however, there exists a mean- 
field on large scales 1$ 100 due to the survey boundary in 
the standard estimator. On the other hand, even without 
source apodization, the bias-reduced estimator suppresses 
the mean-field significantly. If we use the source apodization 
or C _1 filtering, the mean-field is suppressed significantly, 
and the amplitude close to the Monte-Carlo noise level. 

In Fig(3] we show the case with r m = 20'. Even in 
this case, for standard quadratic estimator, either the source 
apodization or C" 1 filtering suppresses the mean field sig- 
nificantly compared to the case without these two filter- 



ing methods, but, similar to r m = 10', there are still large 
mean field at large scales. For the bias-reduced estimator, 
the mean field is suppressed down to the Monte-Carlo noise 
floor. 



4.3 Power spectrum of lensing estimator 

Next, we show the power spectrum of the lensing estimator 
computed from unlensed Gaussian simulations, i.e., the re- 
construction noise bias. The reconstruction noise bias of the 
i-th realization map is computed as 

^'wJ^W* (46) 

and then the mean of 20 realizations is compared with 
the full-sky, diagonal filtering expectation. Note that, for 
the reduced bias estimator, the noise bias is modified, 
the standard reconstruction noise bias divided by 



i.e. 

(1- R*> M R?'*). 

In Fig. |4j we show the reconstruction noise bias, for 
the same cases as shown in Fig. We note that for both 
apodization and C _1 filtering, the reconstruction noise bias 
agrees well with the analytical approximation using either 
the standard estimator (after mean-field subtraction) or the 
bias-reduced estimator. The advantage of the bias-reduced 
approach is that it may be less strongly sensitive on how 
accurately we model the statistics of the underlying fluctu- 
ations. The error bars for the bias-reduced result are also 
smaller than in the standard approach, as this estimator is 
the optimal one in the presence of the statistical anisotropy 
generated by masking. 

In Fig. [5j we show the pseudo-scalar lensing potential. 
The reconstruction noise of the pseudo-scalar lensing po- 
tential is modified and the normalization is biased by the 
mode coupling due to the presence of sky cuts and masks. 
Note that there is no characteristic feature at large scales 
as there is for the scalar-lensing potential. This is because 
the estimator of the pseudo-scalar lensing potential is not 
significantly biased by the masking mean-field, as discussed 
in Sec. [3] Similar to the case with bias-reduced estimator 
of the scalar-lensing potential, using source apodization, the 
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Figure 6. Same as Fig. [4] but using simulated lensed temperature maps. The gray lines show the theoretical reconstruction noise bias, 



N 



0,(0) 



reconstruction noise agrees well with the analytical predic- 
tion for small source holes. 

Finally, in Fig. [6] to see how well the bias can be re- 
duced even in the presence of lensing field, as well as how 
well the estimator normalization is described by the full-sky 
equation, we show the angular power spectrum of lensing 
estimator, computed from Eq. (|46|1 but using the estimator 
reconstructed from i-th realization of simulated lensed map. 
The theoretical prediction is the sum of the reconstruction 
noise, Af (0) =A x e , and the power spectrum of scalar-lensing 
potential, Cf^. The results are similar to that in the case 
with the unlensed Gaussian simulation. 



5 SUMMARY 

We have discussed methods for removing the "mean-field" 
and "reconstruction noise" biases which must be accounted 
for in CMB lens reconstruction. Our approach focuses on 
estimating these biases directly from the data itself, reducing 
our sensitivity to the sky model which is otherwise needed 
to determine them. We performed numerical tests of the 
mean-field reduction approach for several different choices 
of filtering, finding it particularly useful for the reduction of 
the large-scale component of the mean-field. 

In our analysis, we have focused on the temperature 
fluctuations, but upcoming and next-generation CMB ex- 
periments will also provide information on the polarization. 
The polarization anisotropies are a more sensitive probe of 
lensing effects and thus the lensing reconstruction from re- 
alistic polarization maps is also worth investigating. The 
method investigated in this paper may be also applicable 
to the polarization maps, and the usefulness of our method 
to the lensing reconstruction from polarization maps will be 
explored in our future work (Namikawa et al in prep.). 
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APPENDIX A: NUMERICAL SIMULATION OF 
LENSED CMB MAPS 

In this section, we describe our method to generate lensed 
CMB maps, base d on the ray-tracing of large-scale struct ure 
simulations fe.g.. ISato et al. I l2009l ; iTakahashi et al.|[201lh . 

Fig. IA1I shows a schematic picture of our ray-tracing 
simulation. We have multiple lens planes but use a flat-sky 
approximation. The horizontal axis is the comoving distance 
r from the observer. The thick vertical lines are the lens 
and source planes, which are placed at equal-distance in- 
tervals of L, and located at r = (i — 1/2) x L with an 
integer i — 1,2, ••■ ,N. Here, i — N corresponds to the 
source plane. We place the source plane at the last scat- 
tering surface (z s = 1090). The distance to the last scat- 
tering surface is t-lss = r(z s = 1090) = 9900ft _1 Mpc in 
our fiducial cosmological model. We place 19 lens planes 
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Figure Al. A schematic picture of our ray-tracing simulation. 
The horizontal axis is the comoving distance r from the observer. 
The vertical thick lines denote the positions of the lens planes 
and the source plane, which are located at r = (i -1/2) X L with 
i = 1, 2, • • ■ , N. Here, we set N = 20. We have multiple lens planes 
but use the flat-sky approximation. Light rays are emitted from 
the observer and are deflected at the lens planes before reaching 
the source plane. The field of view is 10 X 10 deg 2 . 



up to the last scattering surface, i.e., we set N = 20. 
We determine the interval L from the distance to the last 
scattering surface divided by the number of lens planes, 
L = r L ss/(iV - 1/2) = r-Lss/(20 - 1/2) = 507.6h- 1 Mpc. 
Light rays are emitted from the observer and are deflected 
at each lens plane before reaching the source plane. In our 
simulation, the field of view is 10 x 10 deg 2 . We impose a 
periodic boundary condition on the lens planes. 



Al iV-body Simulations 

In order to obtain the particle distribution and the gravita- 
tional potential on the lens planes, we run iV-body simula- 
tions in a cubic box, and then project the particle positions 
into two dimensions perpendicular to the li ne-of-sight. We 
use t h e numerical sim ulation code Gadget2 (|Springel et al.l 
l200ll ; ISpringel|[2005[ ). We generate the initial conditions 
based on t he second-order Lagrangian perturbation the- 
ory (2LPT; ICrocce et al]|2006l : iNishimichi et al.ll2009l ) with 
the initial linear power spectrum calculated by CAMB. We 
employ 1024 3 dark matter particles in the simulation box 
of L — 507.6 fe _1 Mpc on a side. The initial redshift is 
Zinit = 70, and we dump the outputs (the particle positions) 
at the redshifts corresponding to the positions of the lens 
planes r = L x (i — 1/2), shown in Fig. IA1I The softening 
length is fixed to be 5% of the mean particle separations, 
which correspond to 25 /i _1 kpc. We prepare five indepen- 
dent realizations to reduce the sample variance. 



A2 Ray-tracing Simulations 

We briefly explain the procedure to trace light rays through 
A^-body data and obtain t he maps of the lensing fields on the 
sourc e plane (see also, e.g. JSato et al. 2009; Takahashi et al.l 
2011T) . We use the code RAYTRIX (|Hamana fe Mellieij 



200 If ) which follows the standard multiple lens plane algo- 
rithm. In the standard multiple lens plane algorithm, the 
distance between observer and source galaxies is divided into 
several intervals. In our case, as shown in Fig. lAU we adopt a 
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Figure A2. The power spectrum of the scalar lensing potential 
at the last scattering surface (z = 1090). The dots with error 
bars are our ray-tracing simulation results calculated from the 20 
retaliations of 10 x 10 deg 2 convergence maps. The black curve is 
the theoretical prediction. 



fixed interval whose value is the same as simulation box L on 
a side. Particle positions are projected onto two dimensional 
lens planes (xy, yz ,zx p lanes) every L. Using the Tri angular- 
Shaped Cloud method (jHocknev fc Eastwood! [l988l ). we as- 
sign the particles onto Ng grids in lens planes, then com- 
pute the projected density contrast at each plane. We use 
Ng = 2048 2 throughout this paper. The two-dimensional 
gravitational potential is solved via the Poisson equation 
using a Fast Fourier Transform. Finally, two dimensional 
sky maps of the convergence, shear, and angular positions 
of light rays are obtained by solving the evolution equa- 
tion of Jacobian matrix along the light-ray path which is 
obtained by solving the multiple lens equation. At high red- 
shifts (z > 12), the density fluctuations are very small and 
give only < 10% contribution to the angular power spectrum 
of the lensing potential a t the multipole £ — 100 — 1000 (see 
iLewis fc Challinoi|[2006l ). Hence, at high redshifts (z > 12) 
we do not include the density fluctuations in our calcula- 
tion. We solve the multiple lens equation up to z = 12 and 
further redshifts the light rays are assumed to propagate in 
straight lines. 

We prepare 20 realizations by randomly choosing the 
projecting direction and shifting the two dimensional posi- 
tions. In each realization, we emit 1024 2 light -rays in the 
field-of-view 10 x 10 deg 2 . Then, the angular resolution is 
10deg/1024 ~ 0.6'. 

In order to check the accuracy of our ray-tracing sim- 
ulation, we calculate the power spectrum of scalar lensing 
potential and compare it with the theoretical model. Fig. lA2l 
shows the power spectrum, ^ 4 C^ /2n, as a function of mul- 
tipole £ at the last scattering surface (z = 1090). The dots 
with error bars are the simulation results calculated from 
the 20 realizations of 10 x 10 deg 2 convergence maps. The 
black curve is the theoretical prediction in which the power 
spectrum of scalar lensing potential is given by the pro- 
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jected (three-dimensional) matter power spectrum weighted 
with the radial lensing kernel along the line-of-sight (e.g., 
iBartelmann fc S chneider! 1200 il l. Here, we use the halo- fit 
model (jSmith et alj |2003) to calculate the non-linear power 
spectrum. As clearly seen in the figure, our simulation re- 
sults agree with the theoretical model very well. 



A3 Lensed CMB Temperature Map 

In this subsection, we introduce our procedure for making 
lensed CMB temperature maps. We prepare these maps as 
follows: 

(i) We obtain the power spectrum of the unlensed CMB 
temperature fluctuations using CAMB. 

(ii) We make an unlensed temperature map on a square 
\fAni radian (~ 203 deg) on a side assuming Gaussian fluctu- 
ations based on the unlensed power spectrum. The angular 
resolution of the temperature fluctuations is set to be 10 
deg/1024 ~ 0.6'. We prepare 20 such unlensed maps. 

(iii) Finally, we calculate the deflection angle d(n) at the 
angular position h using the ray-tracing simulation for the 
1024 2 light rays. Then, we obtain the lensed temperature 
map by shifting the positions h — > h + d(n) on the un- 
lensed map, according to Eq. We have 20 lensed CMB 
temperature maps of 10 x 10 deg 2 with 1024 2 grids. 

We calculate the power spectrum from the 20 lensed CMB 
temperature maps, and the result is shown in Fig. IA3I The 
figure shows the angular power spectrum of lensed tempera- 
ture fluctuations as a function of multipole I. The dots with 
error bars are the mean values and the dispersions calculated 
from the 20 realizations. We use so = 0.8 in the apodization. 
The red symbols are the lensed power spectrum, while the 
black symbols are unlensed. The solid curves are the theoret- 
ical prediction calculated by CAMB. Our simulation results 
agree with the theoretical predication very well. 



APPENDIX B: APODIZED NORMALIZATION 

Apodization/masking leads to a change in the angular power 
spectrum of the lensing estimator. Here we derive the ap- 
propriate correction to the normalization of the estimator 
power spectrum. 

Substituting Eq. (|17[) into Eq. ((6J, the estimator in the 
presence of an general window function is given as 



xt = 



d 2 L 



(2tv) 



d 



d 2 £[ 
(2tt) 2 J (2tt) 



(Bl) 



where we define F£ L = Afff iL /2Cf e C*|f ® L] . From the 
above equation, the angular power spectrum of the estima- 
tor becomes 
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Figure A3. The lensed CMB temperature power spectrum. The 
red symbols are the lensed power spectrum, while the black sym- 
bols are unlensed one. The dots with error bars are our simulation 
results calculated from the 20 realizations of 10 X 10 deg 2 lensed 
maps. Here, we use so = 0.8 in the apodization. The solid curves 
are the results from CAMB. 



where we use Q* L — Q-l- Note that the statistical 
anisotropy of lensed fields restrict the integral of Fourier 
modes, £1, £2, £3 and £4,, to the case with £i+£2—£3 — £i = 0, 
the delta function, 8i 1+ i 2 ^i 3 ^i i , is multiplied in the inte- 
grand of the above equation. If the window functions behave 
as the delta function, the above equation reduces to 



<N 2 > 



d 2 L 

(27V) 2 
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X {QL@t-L@-L@-t+L)W4. (B3) 

The quantity W4 is defined by 
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Eq. (|B3|I means that the power spectrum of reconstructed 
estimator from finite-size masked map is equal to that from 
ideal map multiplied by the quantity W4. 
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